# Forelesning 5, onsdag 25.01.2017

# Repetisjon
# Multippel regresjon - intro
install.packages("ISwR")
library(ISwR)
data(cystfibr)
attach(cystfibr)
#Kolinearitet
cor(cystfibr) #Forklaringsvariablene Age, Weight og Height er korrelerte
multregr <- lm(pemax ~ age + weight+height) #Hvorfor er ikke age lenger signifikant?
summary(multregr)

#1# To kategorier kodet ved en dummy variabel - sex
multregr <- lm(pemax ~ age + sex)
summary(multregr)



#4# En enkel balansert design
x1 <- rep(1:5,5)
x2 <- rep(1:5,c(5,5,5,5,5))
cor(x1,x2) 
y <- 0.5*x1 + 0*x2 + rnorm(25,sd=1) # samme model y~x1 som over
balans <- lm(y~x1+x2)
summary(balans)
summary(lm(y~x1)) # betahat_1 forble den samme også etter at x2 ble fjernet. Hvorfor?




#5# Fisher's F-fordeling
par(mfrow=c(2,2))
curve(df(x,df1=46,df2=24),0,5)
curve(pf(x,df1=46,df2=24),0,5)
uu <- rchisq(1000,df=46)  
vv <- rchisq(1000,df=24)
hist((uu/46)/(vv/24),freq=F) #Hva gjør vi her?
par(mfrow=c(1,1))



#6# To utvalg: Like varianser?
heights <- read.table("http://www.math.ntnu.no/~jarlet/statmod/heights.dat")
attach(heights)
var(offspring[sex==0]) # test utført "manuelt"
var(offspring[sex==1])
2*pf(23.26/40.58,46,24) # p-verdi
var.test(offspring[sex==0],offspring[sex==1])
detach(heights)

#7# Enveis variansanalyse
## Eksempel 1
library(ISwR)
data(red.cell.folate)
attach(red.cell.folate)
View(red.cell.folate)
stripchart(folate~ventilation,vert=T)

## Tilpasser modellen og lager en variansanalysetabell
folmodl <- lm(folate~ventilation)
anova(folmodl)
summary(folmodl)  #Hvilken gruppe er forskjellig?

pairwise.t.test(folate,ventilation,p.adj="bonferroni")

## Gjør det samme "for hånd" ...
grandmean <- mean(folate)
groupmeans <- tapply(folate,ventilation,mean)    # gruppe gj.sn., ventialtion kategorisk
ngroup <- tapply(folate,ventilation,length)  # antall obs i hver gruppe

SSD_T <- sum((folate-grandmean)^2)
SSD_B <- sum(ngroup*(groupmeans-grandmean)^2)
SSD_W <- sum((folate - groupmeans[ventilation])^2)
SSD_T
SSD_B
SSD_W
MS_B <- SSD_B/(3-1)
MS_W <- SSD_W/(22-3)
MS_B
MS_W
curve(df(x,df1=2,df=19),0,10)  # plott for F-fordelingen under H_0
abline(v=MS_B/MS_W) # observert F-verdi
detach(red.cell.folate)

#8# Toveis variansanalyse
library(ISwR)
attach(heart.rate)
View(heart.rate)
anova(lm(hr ~ subj+time)) 
sum((hr-mean(hr))^2) # Total sum of squares SSD_T
interaction.plot(subj,time,hr)
interaction.plot(time,subj,hr)  #Hvilket av plottene ser best ut?
detach(heart.rate)

#9# Treveis variansanalyse
heli <- read.csv("http://www.math.ntnu.no/~diserud/ST2304/helicopterdata.csv")

summary(heli)
#heli$flighttime <- rnorm(24,mean=10)   #Kan evt. prøve å erstatte flighttime med denne - hva skjer da i variansanalysen?
attach(heli)
anova(lm(flighttime ~ size + wing + clip))
summary(lm(flighttime ~ size + wing + clip))
anova(lm(flighttime ~ wing)) 
summary(lm(flighttime ~ wing)) # tolkning av dummy koding av faktoren wing
model.matrix(flighttime ~ wing) # dummy variablene som brukes kan ses her...
cbind(heli[,c(2,4)],model.matrix(flighttime ~ wing)[,-1]) # vises ved siden av originaldata
detach(heli)